Multiple time scales hidden in heterogeneous dynamics of glass-forming liquids 
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A multi-time probing of density fluctuations is introduced to investigate hidden time scales of 
heterogeneous dynamics in glass-forming liquids. Molecular dynamics simulations for simple glass- 
forming liquids are performed, and a three-time correlation function is numerically calculated for 
general time intervals. It is demonstrated that the three-time correlation function is sensitive to 
the heterogeneous dynamics and that it reveals couplings of correlated motions over a wide range of 
time scales. Furthermore, the time scale of the heterogeneous dynamics Thetero is determined by the 
change in the second time interval in the three-time correlation function. The present results show 
that the time scale of the heterogeneous dynamics r^tero becomes larger than the a-relaxation time 
at low temperatures and large wavelengths. We also find a dynamical scaling relation between the 
time scale Thetero and the length scale £ of dynamical heterogeneity as Thetero ~ £ z with z — 3. 

PACS numbers: 61.43.Fs, 64.70.P-, 61.20.Lc 



The understanding of drastic slowing down and non- 
exponential relaxations in glasses and dense colloids and 
the jamming of granular materials is one of the most 
challenging problems in condensed matter physics. Ex- 
tensive studies have been carried out through experi- 
ments, computer simulations, and theories [3, 0|. Re- 
cently, the concept of "dynamical heterogeneity" in glass- 
forming liquids has attracted much attention and has 
been considered as an essential notion for understand- 
ing the slow dynamics. Near the glass transition tem- 
perature, the dynamics becomes spatially heterogeneous 
and include the coexistence of mobile and immobile cor- 
related regions. A number of molecular dynamics sim- 
ulationsjajl i, i, 0, 0, i, G3, El El, El and experi- 
ments [l5 1 have detected the existence of dynamical 
heterogeneity with various visualizations. 

To characterize the heterogeneous dynamics in glass- 
forming liquids, explorations that can provide more de- 
tailed information that is not available from conven- 
tional two-point correlation functions, are essential and 
significant. In fact, the growth of the length scale 
£ with decreasing temperature has been measured in 
terms of four-point correlation functions by simulations 
0, 1, H EL El IH EE E3, EI , experiments 0, M S3 , 
and mode-coupling theory [22| . 

Another important issue is the temporal details of the 
dynamical heterogeneity, such as the lifetime and the re- 
location time, i.e., the time scale over which slow (fast) 
moving particles remain slow (fast). The central question 
is whether or not the lifetime of dynamical heterogeneity 
Thetero is comparable to the so-called a-relaxation time r a 
determined by the two-point correlation function [19L |20I| . 
If the deviation between the two time scales becomes 
large near the glass transition temperature, the details 
of the relaxation processes of spatially heterogeneous dy- 
namics are essential to understanding the drastic slowing 
down. 

In order to quantify Thetero, a multiple time exten- 
sion of the density correlation function is necessary, since 
the two-point correlation function averaging over the full 



ensemble cannot distinguish among dynamics of sub- 
ensembles. This idea has been applied to various exper- 
iments such as multidimensional nuclear ma gnet ic reso- 
nance, hole-burning, and photobleaching (l9l |20L l23l . l24l . 
l25l . [2q . Recently, several experiments have provided the 
evidence to support that the lifetime of heterogeneous 
dynamics Thetero becomes larger than the structural re- 
laxation time near the glass transition (2f| SH . In com- 
puter simulations, the multiple time extension has also 
been employed to determine the lifetime of heterogeneous 
dynamics [| ©, ES S3, SI HI M, SH S3 • In these calcu- 
lations, however, due to intense calculations several time 
intervals are fixed at a characteristic time scale, and thus 
only limited information on the lifetime of the heteroge- 
neous dynamics has been provided. 

In this paper, we present the comprehensive informa- 
tion regarding the lifetime of dynamical heterogeneity 
and its temperature dependence. We perform molecu- 
lar dynamics (MD) simulations for simple glass-forming 
liquids at various temperatures and investigate the dy- 
namics in terms of a multi-time correlation function that 
shows the coupling of motions with various time scales 
in the heterogeneous dynamics. The wave vector depen- 
dence of the lifetime of heterogeneous dynamics is also 
studied. Furthermore, the dynamical scaling relation be- 
tween the lifetime and the length scale of the dynamical 
heterogeneity is presented. 

Similar to earlier studies [H SI, SI, S3, S3 , we define 
the three-time, i.e., four-point, correlation function with 
times 0, n, T2, and T3, 
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\ j=i 

(1) 

where 6Fj(k, 0, r) = cos[fc • Arj(Q, t)] — F s (k, r) is the in- 
dividual fluctuation in the incoherent intermediate scat- 
tering function defined as F 8 (k, r) = ((1/-/V) Y^j=i cos (k ' 
Arj(0,r))) with the wave vector k and k = \k\. Here, 
Arj(0, t) = rj(r) — rj(Q) is the displacement vector be- 
tween two times, and r of j-th particle. The (■ ■ • ) repre- 
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FIG. 1: Two-dimensional plots of the three-time correlation 
function F4(k,ti,t2, £3) of component 1 particles at waiting 
time ti = for various temperatures T = 0.473 (a), 0.352 (b), 
0.306 (c), and 0.289 (d). For each temperature, r a is indicated 
by the arrow. The wave vector k is chosen as k = 27r. 

sents the ensemble average over the initial time t = and 
the angular components of the wave vector. The defini- 
tion of the time interval ij is ti = — t,-_i , where tq = 0. 
It is noted that Fi(k, ti,t 2 , £3) expresses the correlations 
between fluctuations in the two-point correlation func- 
tion F s (k,t) between two time intervals, t\ = t± and 
tj, = T3 — r 2 . If the dynamics are homogeneous and if the 
motions between the two intervals t\ and t% are uncorre- 
cted and decoupled, the multi-time correlation function 
Fi{k, ti,t 2 , £3) should become zero. On the other hand, if 
the dynamics are heterogeneous, spatially heterogeneous 
structures governed by the dichotomy between fast- and 
slow- moving regions lead to finite values of F±(k, ti,t 2 ,t 3 ) 
due to correlations of motions between the two intervals 
t\ and £3. Furthermore, the progressive change in the 
waiting time t 2 = T3 — t 2 of F±(k, t\, t 2 , £3) enables us 
to investigate the time scale of the heterogeneous struc- 
ture in glass-forming liquids. In other words, the first 
two-point function 8Fj(k, 0, n) is capable of selecting 
the sub-ensemble of slow (fast) contributions in hetero- 
geneities and the total function 5Fj(k, 0, Ti)SFj(k, t 2 , r 3 ) 
can reveal how long the difference in the dynamics of 
the sub-ensembles remains over the waiting time t 2 [28| . 
It should be noted that multi-time correlation functions 
are exploited in nonlinear multidimensional spectroscopic 
studies of liquids and proteins to understand the detailed 
dynamics, e.g., the transition from inhomogeneous to ho- 
mogeneous dynamics and the couplings of molecular mo- 
tions Him. 

To calculate our three-time correlation function 
Eq. |T]), we have generated long-time trajectories of MD 
simulations for a glass-forming binary mixture composed 
of soft sphere components 1 and 2. Particles interact via 
a soft-core potential v a b(r) = e(a a b/r) 12 , where a a b = 
(fT a +fJb)/2 and a, b S {1, 2}. The interaction is truncated 
at r = 3(Ti. The size and mass ratio are <J\/cr 2 = 1/1.2 
and rn\/m 2 — 1/2, respectively. The units of length, 
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FIG. 2: Two-dimensional plots of the three-time correla- 
tion function F4,(k,t\,t2,tz) of component 1 particles for 
T — 0.289 at several waiting time ti — (a), 0.5T a (b), r a (c), 
and 2r a (d). r a is indicated by the arrow. The wave vector 
k is chosen as k = 2n. 

time, and temperature are <7i, tq = (miaf/e) 1 ^ 2 , and 
e/ks in this paper. Details of the simulations are given 
in previous papers 0, H|. The systems are composed 
of N = Ni + N 2 = 1000 particles with a fixed den- 
sity p = N/L 3 = 0.8 and a composition Nt/N = 0.5. 
The corresponding system linear dimension is L = 10.8. 
Simulations were carried out at various temperatures 
T = 0.772, 0.473, 0.352, 0.306, and 0.289 with a time 
step At = 0.005. Periodic boundary conditions were 
used in all simulations. At each temperature, the sys- 
tem was carefully equilibrated in the canonical condition, 
and then, data were taken in the microcanonical condi- 
tion. We have defined the a-relaxation time r a for each 
temperature by F s (k — 2ir,T a ) = e _1 for component 1 
particles, where k = 2tt corresponds to the wave vector 
of the first peak of the static structure factor. 

In Fig. [TJ we show the two-dimensional plot of 
Fn(k, ti, t 2l £3) of component 1 particles at the zero wait- 
ing time, t 2 = 0, for various temperatures, k is cho- 
sen as k = 27T. It is demonstrated that the intensity of 
F^{k, t 1: t 2 , i 3 ) grows with decreasing T, due to correlated 
motions of heterogeneous dynamics. This indicates that 
particles located in slow (fast) moving regions during the 
first time interval t\ tend to remain slow (fast) during the 
second time interval £3. The line shape of F^(k, t\,t 2 , £3) 
is seen along the diagonal line t\ = £3, and the time at 
which i<4 has a maximum value is approximately given by 
the a-relaxation time r Q . We here remark that the diag- 
onal part t\ — <3 at t 2 = corresponds to the three-time 
correlation function used in earlier studies [13, [H, [2§| 
to judge whether the relaxation type of the dynamics is 
homogeneous or heterogeneous. 

We next investigate how the heterogeneous dynamics 
evolve with waiting time t 2 . Figure [2] shows that the 
waiting time dependence of Fi(k, tt, t 2 , ta) at the lower 
temperature T = 0.289. In addition, in Fig. [3J we plot 
the diagonal line along t\ = £3 of F4 at various t 2 val- 
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FIG. 3: Diagonal part of the three-time correlation function 
F&(k, t, t%,t) of component 1 particles for T — 0.289 at vari- 
ous waiting times t% = 0, 0.5r a , T a , 2r Q , 3r a , 4r Q , 6r Q , and 
10r a from top to bottom. The wave vector k is chosen as 
k — 2tt. The Q-relaxation time and the times at which the 
non-Gaussian parameter <X2(t) and the dynamic susceptibil- 
ity Xi{k = 27T, t) reach their peak values are indicated as T a , 
tngp and r X4 , respectively. 

ues. It is seen that the correlation gradually decays with 
increasing the waiting time ti and tends to become zero 
for t2^oo- The presence of motions correlated with 
the a-relaxation time scale is clearly observed, even for 
t2 = 10r a . Moreover, the off-diagonal part of F\ ex- 
ists as i 2 increases. This indicates that motions between 
the a-relaxation and other relaxations with different time 
scales, such as /3-relaxation, are also coupled, even for 
large ti values. Here, it is seen in Fig. [3] that the peak 
position of F4 changes slightly with ti. For t% = 0, 
the peak appears around the time tngPj where the non- 
Gaussian parameter a 2 (i) = 3(Ar 4 (£))/5(Ar 2 (i)) 2 — 1 
has a peak. On the other hand, for large we find 
that the peak of F4 tends to shift to the slower time r X4 , 
where the non-linear dynamical susceptibility defined as 
Xi(k,t) = N([(l/N) Ejli^KM,*)] 2 ) with k = 2tt has 
a peak. 

Here, we determine the lifetime of heterogeneous dy- 
namics Thctcro and examine its temperature dependence. 
As in the previous studies [U, [HJ , we define the volume 
of Fi(k,ti,t2,t3) by integrating over t\ and t% as 

A hct cro(M 2 ) = / j Fi{k,h,t 2 M)dtxdh. (2) 

Figure H shows A ho tcro(fc, h)/ ^heteto{k, 0) with k = 2tt 
as a function of the waiting time i 2 normalized by the 
a-relaxation time r a . As seen in Fig. [J] Ahctcro rapidly 
decays to zero at higher temperatures, and the time scale 
is comparable to r a . In contrast, at lower temperatures, 
the relaxation of Ahetero occurs on a time scale larger 
than r a . A llc tcro(M2)/A ho tero(fc,0) can be fitted by the 
stretched-exponential function exp[— (£2 /Thetero) c ], where 
"Hictcro can be regarded as the lifetime of the heteroge- 
neous dynamics. We plot Thctcro for each temperature T 
and wave vector k in Fig. \5\ In Fig. it is found that 
Thctcro becomes large as the temperature T decreases. In 
practice, the lifetime Thctcro is about 6t q with c ~ 0.5 at 
T = 0.289 and k = 2tt. Furthermore, we confirm that 
Thetero systematically increases with decreasing fc, imply- 
ing that couplings of large scale motions have long time 




FIG. 4: Waiting time t?, dependence of the integrated three- 
time correlation function Ahotero(fc, i2)/Ah e tero(&, 0) for vari- 
ous temperatures. The waiting times are normalized by the a 
relaxation time r a for each temperature. The wave vector k 
is chosen as k — 2n. The solid curve is determined by fitting 
with the stretched-exponential form for each temperature. 

scales in heterogeneous dynamics. Such strong deviations 
and decouplings between Thctcro and r a with decreasing 
T have been observed in experiments (25l . |26| and simula- 
tions A more recent study has demonstrated 
that the time scale for Fickian diffusion increases faster 
than r a [3(|, which would be related to our results. An- 
other recent study has also mentioned the existence of 
the slower relaxation time is not r a , but the lifetime of 
dynamical heterogeneity [37| . 

Finally, it is of great interest to examine the dynam- 
ical scaling relation between the time and length scales 
of dynamical heterogeneity. As is mentioned, the four- 
point correlation function is used to extract the length 
scale £ of the dynamical heterogeneity, where the fluctu- 
ation in the two-point correlation function can be con- 
sidered as an order parameter as seen in critical phe- 
nomena. Various dynamical scaling relations have been 
proposed 0, H E3] j however the time scale is practically 
chosen as the relaxation time of the two-point correlation 
function, r a . Here, we note that the time scale associated 
with £; should be the average lifetime of the fluctuations 
in the two-point correlation function, i.e., Thctcro- From 
the inset of Fig. we observe a quantitative power law 
behavior between Thctcro and r a as Thctcro ~ with 
C = 1.5 for k = 2tt. A previous study reports that 
a scaling relation between r a and £ as r a ~ £ z with 
z' = 2 in the present binary soft sphere system [7| . Note 
that in Ref. Q the correlation length £ is estimated in 
terms of the static structure factors of bond breakage 
processes among adjacent particle pairs, which is essen- 
tially the same as the £ determined by the four-point 
correlation function consis ting of fluctuations in the lo- 
cal mobility [TTL [T^. IT3L fT^. [171] . Employing this relation, 
we obtain a new dynamical scaling relation as Thetero ~ £ z 
with z = z' Q = 3, which characterizes the relaxation pro- 
cesses of the dynamical heterogeneity. 

In summary, we have investigated a multi-time cor- 
relation function to quantitatively characterize the time 
scale of dynamical heterogeneity. The two-dimensional 
plots reveal couplings of particle motions with various 
time scales, and progressive changes in the waiting time 
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FIG. 5: Lifetime of heterogeneous dynamics Thetcro normal- 
ized by the a-relaxation r Q versus temperature T for ft = ir 
(circles), 2n (squares), 3n (diamonds), and 4ir (triangles). In- 
set: Relation between Thetcro and r a for k — 2tt. The straight 
line with the slope 1.5 is a viewing guide. 

allowed us to extract the time scale of heterogeneous dy- 
namics in glass-forming liquids. It is demonstrated that 
the lifetime of heterogeneous dynamics becomes rapidly 
larger than the a-relaxation time as the temperature de- 
creases, which is compatible with recent optical exper- 



imental results I25l. |26[ and computer simulation stud- 
ies [1, HH, HH, I36l.l37l] . We have also studied the wave vec- 
tor dependence of the lifetime, which is found to system- 
atically increase for large length scales. We have further- 
more confirmed the dynamical scaling relation between 
the time and length scales of dynamical heterogeneity as 
Thetcro ~ £, z with z — 3, where both variables are deter- 
mined in terms of the four-point correlation function. It 
is worth mentioning that the existence of the slower time 
scale Thetero than r a may give the new insights into the 
characterizing of the violation of the Stokes-Einstein re- 
lation P. l35l l38l and non-Newtonian behaviors in sheared 
glassy liquids!?], [39[ , as is discussed in Ref . [37j • 
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